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Abstract. We employ Monte Carlo simulations to numerically study the temporal evolution and transient 
oscillations of the population densities, the associated frequency power spectra, and the spatial correlation 
functions in the (quasi-) steady state in two-dimensional stochastic May-Leonard models of mobile individ- 
uals, allowing for particle exchanges with nearest-neighbors and hopping onto empty sites. We therefore 
consider a class of four-state three-species cyclic predator-prey models whose total particle number is not 
conserved. We demonstrate that quenched disorder in either the reaction or in the mobility rates hardly 
impacts the dynamical evolution, the emergence and structure of spiral patterns, or the mean extinction 
time in this system. We also show that direct particle pair exchange processes promote the formation of 
regular spiral structures. Moreover, upon increasing the rates of mobility, we observe a remarkable change 
in the extinction properties in the May-Leonard system (for small system sizes): (1) As the mobility rate 
exceeds a threshold that separates a species coexistence (quasi-) steady state from an absorbing state, the 
mean extinction time as function of system size N crosses over from a functional form ^ e^^ /N (where 
c is a constant) to a linear dependence; (2) the measured histogram of extinction times displays a corre- 
sponding crossover from an (approximately) exponential to a Gaussian distribution. The latter results are 
found to hold true also when the mobility rates are randomly distributed. 
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have little effect on the dynamical evolution. In this paper, 
our goal is to numerically study the effect of spatial disor- 
der on species coexistence in two-dimensional stochastic 
May-Leonard model variants. To this end, we implement 
the reaction (and mobility) rates for the May-Leonard dy- 
namics as quenched random variables that for each lattice 
site are independently drawn from a truncated Gaussian 
distribution, and subsequently held fixed during the sim- 
ulation run. We here explore (i) the self-organization of 
the population in the coexistence phase, (ii) compute the 
spatio-temporal correlation functions, and (iii) investigate 
the statistics of species extinction times (for small system 
sizes). Our main results can be summarized as follows: 

(1) We demonstrate that quenched spatial disorder 
has only minor effect on species coexistence in the May- 
Leonard model, which together with the results reported 
in Ref. [28 , shows that RPS models (in the presence or 
absence of total particle number conservation) form a class 
of systems that are robust against environmental variabil- 
ity. Remarkably, this statement is true even when spatial 
disorder affects the particles' mobility which is known to 
drastically impact species coexistence. 

(2) We study the combined effect of pair exchange and 
hopping processes, and demonstrate that the former are 
more important for the formation of robust spiral wave 
structures. 

(3) We compute the extinction times, defined as the 
time when the first one of three species dies out, in (small) 
spatially extended systems. We thus find that the mean 
extinction time (MET) increases sharply with system size 
N when the mobility rate is low and the system is in 
the (long-lived met ast able) coexistence state. However, 
once the mobility exceeds the threshold beyond which 
species extinction is prevalent, the MET function switches 
to a linear dependence on N. Correspondingly, the extinc- 
tion time distribution is found to cross over from approx- 
imately exponential with prominent tail at large times, to 
a bell-shaped near-Gaussian function. 

This paper is structured as follows: In Sec.[2l we define 
the stochastic spatial May-Leonard model as a four-state 
spatial rock-paper-scissors (RPS) game without conser- 
vation law for the total particle number, and briefly dis- 
cuss the well-established results from the mean-field rate 
approximation approach. In Sec. [3l we first explain how 
this stochastic spatial system is implemented in our Monte 
Carlo simulations on a two-dimensional lattice, and intro- 
duce the quantities of interest. Then we present and ana- 
lyze our simulation results. Finally, in Sec. H] we conclude 
with a discussion and interpretation of our findings. 



2 Model and rate equations 

In the mean-field approximation, our spatial system re- 
duces to the original May-Leonard modeQ [5 . We let all 
populations live on a square lattice, with each lattice site 



^ At the mean-field level, the model considered here corre- 
sponds to the original May-Leonard model [5] with parameters 
a = 1 and /3 = 1 + cr//i, and time measured in units of If/a. 



occupied with at most a single individual. We therefore al- 
low four states per site: three interacting particle species 
that we label A, 5, and C, and an empty state 0. The 
model is defined through the following set of binary preda- 
tion and offspring production reactions between the three 
particle species [5) [T5) [T6 ] : 



AH 

C- 
X 



B A with rate a ; 

C ^ + 5 with rate a ; 

A ^ + C with rate a ; 

- ^ X + X with rate /i , 



(1) 

(2) 



where X G (A, 5, C) refers to any one of the three species. 
Note that in contrast with the conventional rock-paper- 
scissors model [T8l[2Ql[28l[29l[3Ql [3T1. the total particle num- 
ber is not conserved by these reactions, owing to the sep- 
aration of predation and reproduction processes. In ad- 
dition, in our spatially-extended system we consider the 
nearest-neighbor particle exchange and hopping processes 
on a two-dimensional square lattice (with periodic bound- 
ary conditions; here again X, F = A, 5, C): 

X + r ^ r + X exchange rate e ; (3) 
X + ^ + X hopping rate D . (4) 

It is worth mentioning that the models studied in Refs. [151 
[T6l [T7] do not separate the hopping process (jH from pair 
exchange (jSj); i.e., e = D. Therefore, when letting e = D 
and in the absence of any quenched disorder, our spatial 
model coincides with the one investigated in Refs. [T5j[T6j 
[IT]. 

May and Leonard [5 studied the associated determin- 
istic mean-field rate equations and obtained the temporal 
evolution of the population densities. Let a(t), 6(t), and 
c{t) represent the population densities (concentrations) of 
species A, 5, and C, respectively. Since at most one indi- 
vidual is allowed on each site in the simulation, within the 
mean-field approximation, the overall population density 
p{t) = a{t)-\-b{t)-\-c{t) restricts the reproduction processes 
dU. Therefore, the corresponding rate equations are 

dta{t) = a{t)[p{l-p{t))-ac{t)] , 

dtb{t) = b{t)[p{l-p{t))-cja{t)] , (5) 

dtc{t) = c{t)[p{l-p{t))-ab{t)] . 

The coupled rate equations (j5j) yield four linearly un- 
stable absorbing states (a, 6, c) = (0, 0, 0), (1, 0, 0), (0, 1,0), 
and (0,0,1), and one reactive fixed point (a*,6*,c*) = 
^(1,1,1), where p* = ^^^^ , representing coexistence be- 
tween the three species. Linearizing around the coexis- 
tence fixed point leads to 




= L 



5a 



5c 



(6) 



where 5a{t) = a(t) - a*, 5b(t) = b{t) - 6*, and 5c{t) = 
c{t) — c*, and with the linear stability matrix L 



L 



3p 



p p p^ (J 
p -\- a p p 
p p-^ a p 



(7) 
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Its eigenvalues are Ai = — /i, and A2,3 = 2(3^+cr) [-*- ^ V^^]' 
which demonstrates that the fixed point is locally stable 
only in one direction of parameter space (the eigenvector 
associated with the negative eigenvalue Ai), and generally 
linearly unstable. As elaborated in Refs. [15f il6^ il7j, the 
system dynamics quickly approaches an invariant mani- 
fold associated with the rate equations (|5|). In the neigh- 
borhood of the unstable interior fixed point (a*, 6*, c*), the 
invariant manifold is tangent to the plane normal to the 
eigenvector of L associated with Ai [17 . On this invariant 
manifold, the trajectories approach the absorbing bound- 
aries of the phase portrait where they linger and form a 
heteroclinic cycle [2',"5^. In this case, any chance fluctua- 
tions can cause species extinction by deviating the trajec- 
tories toward the absorbing boundaries. From the imagi- 
nary part of the complex conjugate eigenvalues A2,3, we in- 
fer the characteristic oscillation frequency uj = ap^ I2\fz. 

3 Monte Carlo simulation results for the 
spatially-extended May- Leonard model 

We investigate the two-dimensional May-Leonard model, 
i.e., the four-state stochastic RPS game defined by the 
reactions (|l|2p (which do not conserve the total parti- 
cle number) on a two-dimensional lattice (typically with 
AT = 256 X 256 sites) with periodic boundary conditions, 
subject to the nearest-neighbor exchange (|3]) and hopping 
(jl]) processes. In order to mimic finite local carrying ca- 
pacities, we impose a maximum occupancy number of one 
particle (of either species) per lattice site. When investi- 
gating the effect of spatial disorder on the model, we will 
treat one of the rates (/i, a, e, or D) at each lattice site 
as a random number drawn from a normalized Gaussian 
distribution truncated at one standard deviation on both 
sides. For example, \i ~ N{m^ n) implies that the rate ji is 
picked from a truncated normal distribution on the inter- 
val [m — n^m^n]^ centered at the value m with standard 
deviation n < m. In practice, a value of ji is drawn from 
A^(m, n) for each site on the lattice and attached to the 
corresponding site at the beginning of each single Monte 
Carlo run. The rate values remain unchanged for all sites 
until the next run is initiated. Therefore, in our model, 
the randomized rates pertain to the lattice sites; but are 
identical for any individual landing on a given site for each 
single run. 

At each simulation step, an individual of any species 
on the lattice is selected randomly; then one of its four 
nearest-neighbor sites, which might be empty or occupied 
by one particle of either three species, is selected at ran- 
dom. Subsequently, the particles undergo the pair reaction 
([T]), reproduction (j2j), exchange (j3j), or hopping (jH pro- 
cesses, according to the respective associated rates. Once 
on average each of the P individual particles on the lattice 
has had a chance to react, reproduce, exchange, or move, 
one Monte Carlo step (MCS) is completed; the infinitesi- 
mal simulation time step is thus 5t ~ . 

When studying the effect of quenched spatial disorder 
in the reaction and mobility rates, we shall focus on inves- 



tigating a base model with (average) mobility rate set to 
5, since the corresponding pure system displays clearly es- 
tablished spiral waves (see Fig. Ilbl below). Whenever nei- 
ther fixed rate values nor their distribution are specified 
below, the following default rate values were implemented 
in the simulations: /i = <j = 1, and e = D = 5. We shall 
characterize the emerging spatial structures through in- 
stantaneous snapshots of the particle distribution in the 
lattice, and will depict the temporal evolution of popu- 
lation (spatially averaged) densities. Because of the un- 
derlying symmetry among the species A, 5, and C, one 
representative population suffices and we here report re- 
sults for the spatial average of local population number 
^A{j^t) of species A, i.e., the spatially averaged density 
a{t) = {riAij.t)) = ^ nA{j,t), where j represents the 
site index. We shall also obtain the associated Fourier 
transform 




and compute the equal-time two-point correlation func- 
tions in the (quasi-)steady state, 

CAB{x,t) = {uaU ^ x,t)nB{j,t)) - a{t)b{t) , (9) 

where j denotes the site index, and similarly for Caa{x^ t), 
etc. For our typical system size of = 256 x 256 sites we 
never observed the extinction state in our simulations. In 
fact, as discussed in Ref. [15 and below (see Sec. 3.3), in 
this case one expects the extinction time to grow expo- 
nentially with the system size. Therefore, in order to ac- 
cess absorbing states and numerically measure the mean 
extinction time (MET) Tex = (^ex), where Tex is the ex- 
tinction time for a single Monte Carlo run, and extract 
extinction time distributions, we must consider small sys- 
tems of sizes = 25 to 225. 



3.1 Self-organization in the three-species coexistence 
phase 

In the first row of Fig. 1, we plot typical snapshots of 
the spatial particle distributions at t = 1000 MCS for 
various exchange and diffusion rates, as indicated, while 
jl = a = 1 are held fixed. We observe in Fig. [Tal that all 
three species coexist and a set of entangled spiral patterns 
forms in the system when the mobility rates are compar- 
atively low. Upon increasing the mobility rates, the spiral 
patterns expand, see Fig. [ibl When the spirals' typical 
size ^ at last outgrows the lattice for large particle mobil- 
ity (compare Fig. [TaHTc|) . we anticipate that the system 
essentially acquires the features of its zero-dimensional 
stochastic counterpart (see Sec. 3.3 and Ref. [15 ). In that 
situation, the system evolves towards one of the three ab- 
sorbing states wherein two species become extinct, and 
the surviving species uniformly fills the lattice (uniform 
phase). In finite systems therefore, there exists a thresh- 
old set by the condition £ ^ L = \fN (in a square lattice) 
that separates the absorbing states from species coexis- 
tence. Since the typic al ex tent of the spirals grow dif- 
fusively as £ >/2e, a/2^, this transition should occur 
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(a) e = 0.1,D = 0.1 



(b) e = 5,D = 5 



(c) e = 25, i:> = 25 




(d) cr - Ar(l,0.5) (e) /i- iV(l,0.5) (f) e,D-iV(5,4) 




(g) e = 0,D = 5 



(h) e = 5, D = 



(i) e = 0, D = 25 



Fig. 1: ( Color online.) Snapshots of the spatial particle distribution at t = 1000 MCS for a system with N = 256 x 256 
sites with equal initial densities a(0) = 6(0) = c(0) = 1/3. If not specified otherwise, the corresponding default rate 
values are implemented /i = cr = l, e = i^ = 5. (red/gray: A, yellow/light gray: blue/dark gray: C, black: empty). 



at some critical value Mc of the scaled effective mobility 
M = 2{e,D)/N [T5]. 

In two-species predator-prey systems, quenched spatial 
disorder in the reaction rates can markedly enhance both 
asymptotic particle densities in two-species predator-prey 
systems [32] ; yet this finding is not corroborated in three- 
species RPS models with conserved total population [28 . 
Thus we next explore the effect of spatial variability in 
the reaction as well as in the mobility rates in the spa- 
tial stochastic May-Leonard model where the conserva- 
tion law for the total particle number has been removed, 
and where the mean-field dynamics is not characterized 
by neutrally stable orbits [2l[5l [T8l [2Ql . To this end, we in- 
troduce quenched spatial disorder by treating the rate on 



each site of the lattice as a random variable drawn from 
a truncated Gaussian distribution. As shown in the snap- 
shots in the second row of Fig. 1 (compare with Fig. [Tb|) . 
the presence of spatial clustering can still be observed even 
though small noisy spiral structures dominate the system. 
Thus, the spatial disorder does not markedly affect the 
formation and occurence of spiral patterns. 

The third row of Fig. 1 shows snapshots of the system 
after removing either the two-particle exchange process 
([3]) or pure nearest-neighbor hopping (H]). When exchange 
processes are not allowed, see Figs. p!g| and [TH cluster for- 
mation is still observed, while the spiral waves become 
rather noisy, and it is worth noticing that one additional 
cluster type consisting of only empty sites appears in the 



Qian He et al.: Coexistence in the two-dimensional May-Leonard model with random rates 



5 



system (discernible as black patches in Figs. [Tg| and [Ti]) , 
with measurable consequences on physical observables, as 
will be discussed below. Moreover, in Fig. [Thl where only 
exchange processes are allowed, the spiral pattern bound- 
aries appear quite distinctly sharp, while the empty sites 
are randomly distributed rather than clustered. 

From these results, we infer that the formation of the 
observed spiral patterns is promoted by pair exchange pro- 
cesses. Furthermore, we have also verified that the above 
scenario is not affected when the rates are randomly dis- 
tributed and therefore remains robust against spatial vari- 
ability of the reaction rates. 



3.2 Time evolution and spatio-temporal correlation 
functions 

In order to quantitatively characterize the properties of 
the system, the influence of quenched spatial disorder, and 
the effect of pure particle pair exchange processes on the 
evolution of system, we next depict the temporal evolu- 
tion of population density a{t) = (n^O, ^)), the associ- 
ated Fourier transforms a(/), and the spatial auto- and 
cross-correlation functions Caa and Cab , respectively, in 
the quasi-stationary state (here, at time 1000 MCS) in 
Fig. 2. Since the reaction processes are symmetric with re- 
spect to the three species, the quantities associated with 
species A suffice to extract the relevant information for 
all three species. As shown in Fig.[2al the population den- 
sity decreases swiftly at the beginning of the simulation 
runs: as the particles are initially randomly distributed 
and fill the entire lattice, predation reactions ([1]) domi- 
nate and deplete the particle density at the beginning. 
However, with the emergence of spiral patterns, these re- 
actions can only take place along the domain boundaries 
of distinct species. The system then evolves towards the 
(quasi-) steady state with population densities ^ 0.26 for 
the set of rates chosen here, consistent with the mean- field 
prediction a* = = 0.25. 

In Fig. [2bl we depict the amplitude of the Fourier sig- 
nal of the population density. The peak in |a(/)| yields a 
characteristic oscillation frequency ~ 0.004. This is almost 
a factor ten smaller than the prediction from the mean- 
field approximation, / = uj/27t = V^/IGtt ~ 0.034, indi- 
cating a strong downward renormalization as consequence 
of spatial fluctuations and correlations, similar to the situ- 
ation in the stochastic Lotka-Volterra model [33 ,34 but in 
stark constrast with the conserved spatially extended RPS 
model [28 . The finite width of the frequency peak in the 
Fourier plot indicates that the population oscillations will 
decay and ultimately cease after a finite relaxation time, 
consistent with the damped density fluctuations visible 
in Fig. [2al Therefore, in the coexistence phase, the sys- 
tem's dynamics is consistent with the mean-field descrip- 
tion, a feature that is however caused by the important 
influence of the particles' spatial mobility: With relatively 
low but still effective (average) mobility rates (on average 
e = I) = 5), the dynamics of the system is dominated 
by local interactions (reproduction and predation) along 



the boundaries of local spiral clusters. Thereby, the coex- 
istence state is maintained for a very long time and si- 
multaneously effective mobility mixes the system well, re- 
sulting in a remarkably faithful description of the system 
through the mean- field approximation. Furthermore, the 
(quasi-) stationary auto- and cross-correlation functions in 
Figs. [2c] and [2dl decrease from their extremal values at van- 
ishing distance to zero within about ^ = 20 lattice sites. 

More importantly. Fig. 2 shows the influence of spatial 
disorder on the physical quantities in the May-Leonard 
system. We find that quenched randomness in the rates 
does not noticeably affect the temporal evolution of the 
population densities, the associated Fourier transform sig- 
nals, or the decay lengths of the auto- and cross-correlation 
functions, irrespective of which rate is taken as random 
variable. This result demonstrates that our previous ob- 
servation in the four-state RPS model with conservation 
law, where we found spatial disorder to have only minor 
effects, is also valid for the three-species May-Leonard sys- 
tem. Therefore, we conclude that predator-prey systems 
with cyclic competition appear to be generically robust 
against random spatial variations in the predation, prolif- 
eration, or mobility rates. 

As demonstrated in Fig. 3, the formation of cluster 
patterns (see the third row of Fig. 1) ultimately renders 
the various observables qualitatively similar to those sys- 
tems in which particles can move only via exchange pro- 
cesses. However, in Fig. 3 we also observe that the spatially 
averaged (quasi-) stationary density a* and the peak in the 
Fourier transform |a(/)| vary according to the degree to 
which nearest-neighbor hopping is included in the model. 
In particular, when particle dispersal happens solely through 
hopping processes (jl]), the asymptotic species densities are 
relatively low, see Fig. [3al with e = 0, I) = 5; compare 
with the two plots for other rate choices, and the cor- 
responding snapshots [Tbl [Tg| and [Thl We attribute this 
lower overall species fitness to the previously noted emer- 
gence of sizeable clusters of just empty lattice sites in the 
absence of pair exchange processes (|3]), visible as small 
black patches in Figs. [Tg| and [Til In fact, these voids ef- 
fectively buffer the three species against the predation re- 
actions ([T]), but also diminish the total area that can be 
saturated by either population, which results in an over- 
all density reduction in the (quasi-)steady state. In con- 
trast, when particle pair exchange processes are included, 
the influence of empty-site clusters is diminished and con- 
sequently the (quasi-)stationary population densities en- 
hanced (see Fig. [3alfor e ^ 0). Furthermore, we observe 
that the characteristic density Fourier peak frequencies 
in systems where nearest-neighbor hopping is allowed are 
renormalized to even lower values than in runs with just 
pair exchange processes, which tend to better mix the sys- 
tem, see Fig. [3bl (compare the characteristic frequency for 
the runs with D = 5 to those with D = 0).ln Figs. [3cl and 
[3dl we observe the presence of low-amplitude population 
(damped) oscillations during the decay of the correlation 
function, which originate from the spiral structures dis- 
played by the system in the coexistence state (under "effi- 
cient" mobility rates) [TT. In the insets of Figs. [3cl and [3dl 
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-a~N(1, 0.5) 
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(c) Correlation function Caa{x). 
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^^^~N(1,0.5) 






— V— 8~N(5, 4) 






D~N(5, 4) 









40 ^ 60 



(d) Cross-correlation Cab{x). 



Fig. 2: (Color online.) Quantitative observables for a stochastic May-Leonard system with N = 256 x 256 sites, 
starting with equal initial densities a(0) = b{0) = c(0) = 1/3 and in the presence of spatial disorder, averaged over 50 



simulation runs. If unspecified, the default rate values /i 
(c) and (d) were measured at t = 1000 MCS. 



a = l^e = D = b were used. The correlation functions in 



we notice that the correlation functions of the model vari- 
ants with only hopping processes (e = 0) decay to zero in 
a less oscillatory manner than the correlation functions for 
the variants that include particle exchange (e 7^ 0). This is 
a consequence of the observed absence of well-defined spi- 
ral structures (see Figs. [Tg| and [Ti]) when the pair-exchange 
rate is too low to efficiently stir the system. 

In summary, the pair-exchange processes suppress the 
presence of empty-site clusters, and despite rendering the 
spiral structures more diffuse (i.e., more "entangled"), they 
have an overall stabilizing effect on emerging spatial pat- 
terns, and consequently promote the fitness and coexis- 
tence of all three subpopulations A, B and C. 



3.3 Mean extinction times and their distribution 

As we have discussed above, the system is efficiently stirred 
when the pair exchange rate e is high enough (independent 



of the actual value of D). In this case, the system is char- 
acterized by a long-lived coexistence state. However, when 
the pair exchange rate is below some critical value, it has 
been shown that the system settles in an absorbing state 
after an observable amount of time p3^. We here revisit 
and extend the analysis of such a scenario that holds true 
for any values of e or D by computing the mean extinction 
time (MET) and the distribution of extinction times. Here, 
for the sake of simpHcity (and without loss of generality) 
we assume that e = D and the mobility rate therefore 
is M = 2e/N . In this setting, our May-Leonard model 
coincides with the variant considered in Refs. [r5j[T6l[T7], 
where it was shown that the critical mobility threshold is 
Mc ~ 4.5 X 10""^ (when /i = a = 1). Indeed, when the ef- 
fective mobility rate M approaches Mc from below, there 
is a cross-over from a coexistence (quasi-) steady state to 
an absorbing state. 

The MET has been computed as the time when the 
first of the three species dies out. Figure 4 illustrates how 
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Fig. 3: ( Color online.) Quantitative observables for a stochastic May-Leonard system with N = 256 x 256 sites, reaction 
rates /i = a = 1, and starting with equal initial densities a(0) = 6(0) = c(0) = 1/3, with different combinations of 
nearest-neighbor particle exchange and hopping processes, averaged over 50 simulation runs. The correlation functions 
in (c) and (d) were measured at t = 1000 MCS. 



this MET exhibits markedly different behavior when the 
effective mobility rate M = 2e/N is swept through Mc. 
As shown in Fig. l4bl when the mobility is relatively weak 
(M = 10~^), the MET increases with the system size 
approximately according to the (zero-dimensional) func- 
tional form Tq^{N) ~ e^^/A/" (where c is a constant) [3F , 
especially for comparatively large values N ^ 200 . . . 600. 
(For smaller systems with N < 200, the data do not fit 
this functional dependence very well, but may also not be 
as statistically reliable.) Yet the curvature of the graphs 
in Fig. Hal decreases upon raising the effective mobility, 
and the functional dependence on system size becomes re- 
placed with a linear form fe^{N) N for M > Mc. That 
is, when the mobility rate is low, the system is dominated 
by local interactions and species extinction is a rare event 
driven by a large fluctuation after an enormous amount of 
time. In this case, the coexistence of the three species cor- 
responds to a metastable state. Interestingly, Fig.Hclshows 
that spatial disorder in the mobility rate M does not qual- 
itatively affect the behavior of the MET: this very same 



scenario applies even when M is randomly distributed. 
This observation further supports the conclusion that spa- 
tial variability in the mobility rates has little effect on the 
dynamical evolution of the system. 

When the mobility rate increases and exceeds the thresh- 
old, the system is regularly driven towards extinction and 
biodiversity is lost, just as predicted by the zero-dimensional 
formulation of the model. Furthermore, the histograms of 
extinction times plotted in Figs. l4dl and Hel show that the 
extinction time distributions (obtained for small systems 
with A" = 20 X 20 sites) correspondingly evolve gradually 
from an approximately exponential (or Poisson) shape, al- 
beit with fat tails, towards a (roughly) Gaussian distribu- 
tion centered at Tex- In addition, the distributions remain 
unchanged even if spatial disorder is incorporated in the 
model (see the insets of Figs. l4dl and He|) . The approxi- 
mate quasi-exponential (or quasi-Poisson) distribution of 
Fig.|4d]is a characteristic feature of systems where extreme 
events occur only after a very long time, and which are 
hence driven by large rare fluctuations [36 . Here, the rare 
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(C) Tex VS. N. 
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(d) e = T> = 0.1; inset: e = T> - N(0.1, 0.05). 




(e) e = T> = 25; inset: e = T> - N{2b, 10). The dashed 
(red/gray) curve depicts a Gaussian fit. 



Fig. 4: (Color online.) (a), (b), (c) Mean extinction time (MET) Tex as function of lattice size A/", for different values 
of the effective mobility M = 2e/N (here, e = D), obtained from averages over 10000 Monte Carlo runs, starting 
with equal initial densities a(0) = 6(0) = c(0) = 1/3, and reaction rates /i = a = 1. The lattice sizes are = 5 x 5, 
7 X 7, 10 X 10, 12 X 12, 15 x 15, 17 x 17, 20 x 20, 22 x 22, 25 x 25 sites, (d) and (e): Histogram of extinction times 
estimated from 10000 Monte Carlo runs, based on a system with lattice size = 20 x 20. The insets correspond to 
the histograms obtained for random mobility rates: e = D ^ A^(0.1, 0.05) in (d) and e — D ^ A^(25, 10) in (e). 



extreme event is extinction of species that were previously 
coexisting in a metastable state for a long time period. On 
the other hand, the approximately Gaussian distribution 
of Fig. [lilis typical of systems where random fluctuations 
are of weak intensity [37] . In our competing three-species 
system, it is associated with the absorbing state that cor- 
responds to species extinction happening within an "ob- 
servable" time range, as in the zero-dimensional counter- 
part of the model ^15[[T6l [T7] . These numerical observations 
support the method suggested in Ref. [15] for identifying 
the nature of (quasi-) steady states in the simulation of 
RPS models: A reasonable criterion is that the system re- 
mains at the steady state if three species still coexist after 
simulation time t ^ otherwise, the system evolves to 
an absorbing state. 



4 Conclusion 

In this paper, we have demonstrated that quenched spatial 
disorder in either the reaction or the mobility rates does 
not significantly affect the temporal evolution, Fourier sig- 
nals, spatial correlation functions, or mean extinction times 
in stochastic spatial May-Leonard models (i.e., four-state 
RPS models without total particle number conservation) 
in two dimensions. In combination with our previous re- 
sults for conserved three-species RPS systems [28j, we con- 
clude that such cyclic predator-prey systems appear to be 
generically robust against spatial variability of the rates. 
Here, the randomized reaction rates remain attached to 
the lattice sites, mimicking environmental variations that 
do not change over time. As there exist a number of sys- 
tems, e.g. in ecology [12] and microbiology [13], where 
the competition among species is cyclic, an important im- 
plication of our findings is that the environmental vari- 
ability of the parameters can essentially be neglected in 
the mathematical description of those systems. In addi- 
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tion, through removing the hopping process by letting 
D = 0, we observe that particle pair exchange processes 
promote the formation of sharp spiral patterns. In our 
spatial stochastic system, we measure the population os- 
cillation frequency to be much lower than predicted by the 
mean- field rate equations, similar to the situation in two- 
species Lotka-Volterra models [33,34 , but in stark con- 
trast with our numerical results for conserved RPS model 
variants [28]. This downward frequency renormalization 
is enhanced by the presence of nearest-neighbor hopping 
processes. Moreover, we find a remarkable gradual trans- 
formation in the dependence of the mean extinction time 
on system size, and the shape of the associated extinc- 
tion time distribution, when the effective mobility rate 
crosses the critical threshold separating the coexistence 
from the absorbing state: When the mobility rate is low, 
the distribution of extinction times is approximately ex- 
ponential, and species coexistence corresponds to a long- 
lived met ast able state. In this case extinction is driven 
by large, rare fluctuations and the mean extinction time 
essentially grows exponentially with the population size. 
Above the critical mobility threshold, the extinction times 
are approximately distributed according to a Gaussian. In 
this situation, the noise is of weak intensity and the mean 
extinction time grows linearly with the population size. 
Interestingly, we find that these results remain valid for 
both non-random as well as for randomly distributed mo- 
bility rates. 

This work is in part supported by Virginia Tech's Insti- 
tute for Critical Technology and Applied Science (ICTAS) 
through a Doctoral Scholarship. We gratefully acknowl- 
edge inspiring discussions with George Daquila, Uli Do- 
bramysl, Michel Pleimling, Matt Raum, and Royce Zia. 
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